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Abstract 

I develop the formalism that allows calculation of beam envelopes 
through a linear accelerator given its on-axis electric field. Space charge 
can naturally be added using SachererJTj] formalism. A complicating 
feature is that the reference particle’s energy-time coordinates are not 
known a priori. Since first order matrix formalism applies to deviations 
from the reference particle, this means the reference particle’s time and 
energy must be calculated simultaneously with the beam envelope and 
transfer matrix. The code TRANSOPTRlEH is used to track envelopes for 
general elements whose infinitesimal transfer matrices are known, and in 
the presence of space charge. Incorporation of the linac algorithm into 
TRANSOPTR is described, and some examples given. 


1 Introduction 


For beamline design, matching, etc., even in the case of strong space charge 
forces, it is often sufficient to calculate the evolution of size parameters such as 
the rms width and bunch length. The full set of such parameters is the 
6-dimensional “a-matrix”; x, P x ,y, P y , z, P z . Conventionally, the momenta are 
divided by the reference particle’s momentum, making them x', y'. AP/T[] 
Parameter 5 is technically time, but converted to distance (bunch length) by 

'Throughout this note, primes denote d/ds 
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multiplying by the speed f3 0 c of the reference particle. Similarly, parameter 6 is 
A E, converted to P z by dividing by (3 0 c. The a-matrix has 21 parameters, as it is 
symmetric: 6 rms sizes on the diagonal, 15 correlation parameters. 

The mathematical fomalism for this technique, including space charge, was 
established by Frank SachererllTII. The space charge forces depend crucially upon 
the bunch dimensions in configuration space, so it is important that these be 
tracked. In other words, it is insufficient to use a formalism that first integrates 
the equations of motion to derive the transfer matrices, and apply space charge 
effects afterwards. Some implementations such as TRANSPORT and TRACE3D 
divide standard elements into (hopefully sufficiently short) segments interleaved 
with space charge “lenses”. This is crude, approximate and non-adaptive. 

TRANSOPTRlUl uses the envelope formalism, but did not until now include the 
case of beam high intensity bunches being accelerated with a RF accelerator. 

This case is of particular importance for modeling the elecron linear accelerator. 
Short accelerator gaps can often be sufficiently well modelled as infinitesimal i.e. 
as a thin longitudinal lens with also a transverse thin lens focusing component. 
But this is insufficient for extended RF devices and as well misses time-of-flight 
effects when the change in velocity of the reference particle is significant. 
Extended DC longitudinal fields have already been added to TRANSOPTR|f3ll; the 
present study aims at adding the AC case. In the DC case, the reference particle’s 
longitudinal coordinates are inferred from the local potential. Not so the AC 
case; the time and energy versus longitudinal location must be tracked in a 
separate calculation and then the equations of first order deviations applied to 
find the transfer matrices and local envelope optics. 


2 Theory 


This follows closely from the previous note on DC axially-symmetric 
longitudinal fields 01. 
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2.1 Hamiltonian 


With the distance along the reference trajectory s as the independent variable, the 
Hamiltonian is 


H(x, P x , y, P y , t, E\ s) 



E - q$ 

c 


2 


77l 2 C 2 


(p x - qA x ) 2 


(Py ~ qAyf 
( 1 ) 


The case of RF axially-symmetric electric field can be handled entirely with no 
electric potential (<f> = 0), and time-varying vector potential. This has been 
presented a number of times in the past (e.g. Chambers|j4l), but we are interested 
in the following more experimentally-useful case: The electric field along the 
axis £(s) has been measured and is therefore known, and the geometry is exactly 
axially symmetric. Rob Ryne 0 has treated this case, and we use his vector 
potential A(x, y, s, t) directly. 


A x 




£'{s) sin(cot + 9) 

2 uj 

£'(s) sin(o;f + 6) 


2 


UJ 
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£"(*) + 



sin (cut + 9 ) 

CO 


( 2 ) 

( 3 ) 

( 4 ) 


This is Coulomb/Lorenz gauge, satisfies Maxwell equations to second order in 
transverse coordinates, gives correct on-axis £ = —dA/dt = £ cos(cuf + 9). 


A priori, we do not know the reference particle’s energy and time coordinates. 
We need these in order to expand about them. They can be found from the 
equations of motion for x = y = P x = P y = 0: 

d E 9H „ . 

—— = = q£ cosfc ot + 9) 

as at 


dt 

ds 


dH 

~dE 


E 

P 


1 

Po c 


( 5 ) 

( 6 ) 


(From here on, I drop the 0 subscript: f3 and 7 are implicitly assumed to be the 
relativistic parameters of the reference particle.) 


These are solved first and give the functions E 0 (s) and t 0 (s) about which t and E 
are expanded: E = E 0 + A E, t = t 0 + At. So we transform the canonical 
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variables t and — E to (At, — A E), using as generating function 


G = 




( 7 ) 


(Check: ^ = —E, d ( d< ^ E ) = At.) The Hamiltonian gets the added terms 

dG _AE + E 0 (s) 
ds /3(s)c 

Then expanding the square root, we get: 


- A tE' 0 (s). 


H At = ( y c ~ Poj ^qA s -AtE' 0 (s)- 


(A E) 2 {P x - qA x ) 2 + (Py - qA y f 

2/? 3 7 3 mc 3 2 P 

( 8 ) 


In expanding P x — qA x , P y — qA y , the time dependence disappears because it is 
higher order: 


{Px - qA x ) 2 


Px - qS- 


! sin(cut 0 + 0) 


xP x + 


u 


( q£' sin(o;t 0 + 9) 


2 


X 


( 9 ) 


and similary for y. The term linear in At in the expansion of A s about t 0 cancels 
the —AtE' 0 (s) term, as it should but there is a remaining term quadratic in At, 
the bunching effect. This leaves 


-qA s -AtE' 0 (s) 


sin(o;fo + 6 ) ( w 2 (At) 2 

<?£- 1 -^- 

u> \ 2 


r 2 q ( p „ Anjujtp + 9) 

4 ( ' c 2 J cu 

( 10 ) 


Notice the first term here and the first term in eqn.[8]depend only on the 
independent variable and not on the 6 dependent ones. Thus these do not affect 
the equations of motion and we ignore them. We have: 


= -f + -APE 

2 Zp 6 ^ 6 mc 6 


2 

r q 


P 2 rP 

J_ 5 - _ nP’T _- 4 - 

+ 2 P 9 2P + 

+ 2 P 9 2P + 
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We defined here T(s) = sin[o;fo( s ) + 0\/ui to clean up the notation a bit. 

Finally, we wish to transform from (At, —A E) to (z, P z ) = (—/ 3cAt , A E/(/3c)). 
(The reason for the sign change is as follows: an early arrival implies At < 0, but 
this means the particle is ahead so z > 0.) The generating function is 


G = -ficAtP z 


( 12 ) 


dG 


(Check: aM 


-A E, tJ'Tt = z.) The term to be added to the Hamiltonian is 

’ o(P z ) ' 


= ±_ 

ds (3 2 (3 2r y 3 


zP z = 


q£C 

/3cPy 


zP z , 


where C = cos(c ot Q + 6). 
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2.2 Hamiltonian 2 


Ryne[0 has a transformation that gets rid of the second derivative of the on-axis 
elecric field. It’s complicated. At the same time he transforms away the adiabatic 
damping; it’s a neat and didactic trick but not strictly necessary for 
computational purposes. It is simple to just use P x , y ,z directly and then just 
rescale by final P at the end. 

But there’s an easy way to get rid of the second derivative: it turns out that the 
vector potential can be simplified if we use a different Gauge. 


I propose the following function 

^(x,y,s,t) 


£' sin (out + 9) x 2 + y 2 
2 ~ u 2 


(14) 
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Add the gradient of this function to the previous vector potential ( 2|3|4 ). This 
zeroes both A x and A y , leaving 


A a = -£{s) 1 - 


uj 2 x 2 + y 2 \ sin (out + 9) 


c 2 4 J oj 

This is considerably simpler, but now there is a scalar potential: 

Q'ty x 2 + v 2 

<3> = — = S' cos(c ot + 6) ——— 

Now if we expand the Hamiltonian, we get a different result: 


p 2 P 2 n 

1 x i y I " 


p j/3\ r 2 


S ‘ = tp + £ + Wc\ £ ' C - £S A) 2 + 

P 2 2 2^^^ zP z q£uS z 2 


2y 2 P /3c 2 7 2 P /3 2 c 2 2 

(C = cos(cufo('S) + 9), S = sin(o;fo( s ) + 8)) 


( 15 ) 


(16) 


(17) 


This is not only much simpler than eqn.[l3](P r and P y have their usual 
definitions, no transverse cross terms, no £"), but has nice intuitive explanations 
for the individual terms. (1) The factor in parentheses is precisely the integrand 
of eqn. 4 of a short note I wrote in 1985® explaining the derivation of the focal 
power of an RF gap, e.g. a buncher. (2) Taking the limit as uj —> 0 reproduces 
precisely the Hamiltonian of the DC accelerator I derived in 2010OJ. Note that 
in that case, S' = —0". 


2.3 Infinitesimal Transfer Matrix 

A convenient and useful way of representing the equations of motion through the 
optical element is the so-called infinitesimal transfer matrix approach®. The 
infinitesimal transfer matrix F(s) is defined as (T — I)/ds where T is the 
transfer matrix from s to s + ds and / is the identity matrix. With this definition 
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one has for individual particles 


dx 

ds 


= Fx, 


where 


x = 


( X \ 

p x 
y 
Py 

z 

V P, ) 


(18) 


Beams of particles are conveniently represented by the so-called a-matrix; the 
elements of which represent second order moments of the beamflTI. The a -matrix 
and the transfer matrix M transform through the system according to the 
equations 


da 

ds 

dM 

ds 

where F T is the transpose of F. 


Fa - aF T , 

(19) 

FM. 

( 20 ) 


Now that the Hamiltonian for linear motion (eqn.|T7|) has been obtained, it is a 
simple matter to find the infinitesimal transfer matrix F. Writing the equations of 
motion (x' = dH/dP x , P' x = —dH/dx, etc.) in the form of Eqn. 18 the 
following F-matrix is found for the axially symmetric linear accelerator: 


F 


1 0 T o o 0 

.A(s) 0 0 0 0 

0 0 0 T o 

0 0 .A(s) 0 0 

0 0 0 0 | 

v 0 0 0 0 B(s) 


0 N 
0 
0 
0 


i 

7 2 P 



where we have defined: 

As) 

B(s) 


-4- (s'C-SS 

2 Pc \ 
q£uS 
f3 2 c 2 ' 



) 


( 21 ) 


( 22 ) 

(23) 
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2.4 Space Charge 


Space charge forces depend recursively upon the a-matrix elements, and are 
simply added to the focusing elements F 2 n,2m-i\n,m=i,2,3 of the element’s 
infinitesimal transfer matrix such as eqn.| 2 T| above. This technique is used in the 
code TRANSOP TR, as described by de Jong[[2j]. The resulting equations can only 
be solved numerically. 

The given references fll] 0 treat space charge in the non-relativistic regime, so it 
was not obvious that TRANSOPTR was correct in the relativistic regime. There 
are two effects that need to be considered to generalize the equations: the space 
charge magnetic field, and bunch length contraction. For detailed derivations, the 
interested reader is referred to the Ph.D. thesis of FubianiQ. The first effect 
requires dividing the space charge force by y 2 . The second requires that the 
Carlson elliptic integrals’ arguments be modified. The 3 integrals for the 3 major 
axes are 

3 f°° dt 

R D \U, v,w) — - / - 

^ “'° (t + w) yj(t + u) (t + v) (t + w) 

where ( u , v , w ) are (cr 33 , a 55 , a n ) for the tc-axis, (cr 55 , a n , a 33 ) for the y- axis, 
(cTn, cr 33 , cr 55) for the z-axisj^To be relativistically correct, <755 must be replaced 
by 7V55. (See FubianiQ, Appendix J.) 

An interesting limit is the long bunch, since this can be approximated as a 
continuous beam with current I. First of all, it is clear that for this limit to apply, 
bunch length transverse size is not a necessary condition. Rather, 7 times 
bunch length transverse size. This means that for example a 1 mm long by 
1 mm wide electron bunch is already well into the long-bunch regime with 
energy of 10 MeV. 

Secondly, in the long bunch regime, the Carlson integrals governing transverse 
space charge are oc or the inverse bunch length augmented by the 

factor 7. For this reason, the constant governing space charge force in 
TRANSOPTR in the unbunched case is divided by an extra factor of 7. 

2 Note that this assumes the bunch axes are aligned with the beam direction (s) and the chosen 
transverse orientation. If this is not the case, a rotation must be applied. 
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3 Implementation into T RAN S 0 P T R 


In TRANSOPTR, the momentum is dimensionless and in absence of acceleration 
corresponds to angles (x', ?/, z'). When acceleration occurs, angles are not scaled 
canonical momenta. The simplest implementation is to scale (P x , P y . P z ) by 
initial total momentum P t . Then after integration is complete, the momenta can 
be converted back into angles by multiplying by Pi/Pf. In this case, F (eqn. 21) 
becomes: 


F 


( 0 f 0 0 0 0 N 

^ 0 0 0 0 0 

0 0 0 f 0 0 

0 0^00 0 

0 0 0 0 f ^ 

V 0 0 0 0 -f 


(24) 


In order to find the a-matrix at any point s of a beamline, it is sufficient to 
calculate the transfer matrix to that point. This technique is not useful when 
space charge forces are non-negligible since then the optics themselves depend 
upon the a-matrix. In the past space charge has been added to legacy codes such 
as TRANSPORT by inserting a sufficient quantity of defocussing lenses whose 
strength depends upon the local beam size. This is inefficient at best: it is the 
crudest sort of integration. The best way to do this is as implemented in 
TRANSOPTR: the exact Sacherer 6D envelope equations (19) are integrated with 
a higher order integrator such as the Runge-Kutta. 


In many situations, one needs only the cr-matrix and so it is sufficient to solve the 
36 equations of eqns. 19| 3 However, in general one would like to have the ability 
to fit certain optical characteristics, such as point-to-point imaging. It is therefore 
often useful to also know the transfer matrix to any point. This is in fact the 
incoherent transfer matrix. So in order to cover all cases, TRANSOPTR solves 


3 In fact, since this matrix is symmetric, one can get away with solving only 21 equations. One 
can whittle this down further for uncoupled cases; in fact 4 first order or 2 second order equations 
are often all that is needed for a straight beamline with no coupling and no dispersion. These are 
of course the Kapchinsky-Vladimirsky equations. 
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both eqns.[19]and|20j This comprises 72 equations. Internally, the code uses a 
single 12 by 6 matrix. 

Without time-varying fields, the total momentum of the reference particle is 
always known at any s. In a linac, this is no longer true. Thus, the 
implementation needs to track not only the transfer and a matrix elements, but 
also time and energy. These can be thought of as the 5th and 6th “zero-order” 
coordinates, with respect to which the fifth and sixth phase space coordinates are 
measured. Analogously, in the most general cases, the first through fourth 
“zero-order” coordinates may also not be known. This could happen for example 
if given only the field map of a dipole magnet. Thus there could be 6 
“zero-order” coordinates as a function of s, and the usual 6 phase space 
coordinates are differentials with respect to these. 

I therefore have added a row to the 72-element matrix, making it dimensioned 
13-by-6, and there are 78 ODEs. This would allow at a later date adding 
elements whose reference trajectory is not known a priori. 


3.1 Code 


Here follows the f77 subroutine that takes the 13-by-6 matrix SX and calculates 
the 13-by-6 matrix DSX, which are the derivatives of SX with respect to Z =s. 

Note that SX (13,5) is not time t directly but is ct, and SX (13, 6) is not 
energy directly but is 7 — 1. The F matrix is in fact F directly. 

SUBROUTINE SCLINAC(Z,SX,DSX) 

COMMON/PRINT/IPRINT,IQ1,JQ1,IQ2,JQ2,IQ3,JQ3,IQ4,JQ4 
COMMON/VELEM/VELM(21,9999),NELM 
DIMENSION SX(13,6),DSX(13,6) 

DIMENSION TEMP(6,6 ),ROT3(3,3),ROT6(6,6),SIG(3,3) 
DIMENSION WRK(3) ,EVAL(3) ,RF(6 ,6) ,FT ( 6,6) 

LOGICAL STRATE 
common/splind/klob,kloe 
COMMON/FMATRX/F(6,6) 

COMMON/SCPARM /KSC,ISC 
COMMON/CALLS/NSC 
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COMMON/MOM/P,BRHO,pMASS,ENERGK,GSQ,ENERGKi,charge,bncharge 
COMMON/axezrf/omoc,phase,etot, 

$ ZTBLE(9999) ,EZTBL (9999) ,SPCFE(9 999) ,zse,iez,iaxezrf 

COMMON/axez2/pratio 
REAL KSC 
do i=l,4 
dsx(13,i)=0. 
enddo 
NSC=NSC+1 
pratio=l. 

gamml=sx(13,6) ! gamml=gamma-l 

energk=gamml*PMASS 

if (energk.le.0.)then 

write(6,*)'*********BEAM DECELERATED TO A STOP' 
write ( 6, *)' Check parameter values or limits.'' 
call errout 
endif 

c call spline routine to find the on-axis electric field and derivativ 
call spling(kloe,ZTBLE,EZTBL,SPCFE,IEZ,z-zse,EZ,EZp,EZpp) 
ccc=cos(sx(13,5)*omoc+phase) 
sss=sin(sx(13,5)*omoc+phase) 
c omoc is omega/c 

dsx(13,6)=EZ*ccc/PMASS/100. ! denergy/ds ; EZ is in MeV/m, a 

gamma=l.+gamml 
ETA=SQRT(gamml*(gamml + 2.) ) 

BRHO=ETA*PMASS/CHARGE*3.3356397 
BRHO=BRHO*.001 

P=ETA*PMASS !this P is actually Pc 
gsq=GAMMA* *2 
BETA=eta/gamma 

dsx(13,5)=1./beta ! dctime/ds=l./beta 
gammlold=etot/pmass 

ETAold=SQRT(gammlold*(gammlold+2.)) 

Pold=ETAold*PMASS 
pratio=etaold/eta 

c re-calculate space charge parameter as energy changes, 
c (898755.2 is 1/(4pi epsilon_0) ~ c"2/l.e7 in some units) 
KSC=BNCHARGE*CHARGE* 8 9 8 7 55.2/(ETA**2*PMASS) /pratio 


11 



F(1,2)=pratio 
F(3,4)=pratio 

F(2,1)=(EZ*sss*omoc-EZP*ccc/beta)/(200.*Pold) 
F (4,3) =F (2,1) 
c bpob = beta'/beta 

bpob=dsx(13,6)/(gamma*eta**2) 

F (5,5)=bpob 
F(5,6)=pratio/gsq 

F(6,5)=EZ*sss*omoc/(100.*Pold*beta* *2) 

F ( 6, 6) =-bpob 
IF (KSC.NE.O.)THEN 


0 

c START SPACE CHARGE SECTION: 
STRATE=.TRUE. 

c SIG is real-space alone (3D) 
DO 5 1=1,2 
DO 5 J=l, 2 

5 SIG(I,J)=SX(2*1-1,2*J-1) 


gamma=sqrt(gsq) 


C QUICK AND DIRTY Lorent 

z transform. 

SIG 

(1,3) 

=GAMMA*SX( 

1,5) 

SIG 

(2,3) 

=GAMMA*SX( 

3,5) 

SIG 

(3,1) 

=SIG (1,3) 


SIG 

(3,2) 

=SIG(2,3) 


SIG 

(3,3) 

=GSQ *SX( 

5,5) 


c write(36 ,*)' gsq',gsq,sig (3,3) 


c 'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 
c Test whether bunch axes are along coordinate axes 

C ****i<i<*i<***'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

IF (SIG(1,2).NE.0..OR. 

, SIG(1,3).NE.0..OR. 

, SIG(2,3).NE.0.)THEN 

CALL PRESET(ROT6,36,0.) 

STRATE=.FALSE. 


c Find rotation required 

call jacobi(sig,3,3,eval,rot3,nrot) 
c Fill 6D rotation matrix 
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DO 55 1=1,3 
DO 55 J=1,3 

ROT6(2*1-1,2*J-l)=ROT3 (I,J) 
55 ROT6(2*1,2*J)=ROT3(I,J) 

ELSE 

c Bunch is along axes 
DO 52 1=1,3 

52 EVAL(I)=SIG(I,I) 

END IF 


do i=l,3 

if (eval(i).le.1.e-6)then 
eval(i)=1.e-6 
endif 
enddo 

c RD is the Carlson elliptic integral R_D. 

c Scale the arguments. Carlson routine is better if args are order 1. 
scalrd=eval(1) 

if(eval(2).gt.scalrd)scalrd=eval(2) 
if(eval(3).gt.scalrd)scalrd=eval(3) 
rscalrd=scalrd* *(-1.5) 

ELl=rscalrd*RD(EVAL(2)/scalrd,EVAL(3)/scalrd,EVAL(1)/scalrd) 
EL2=rscalrd*RD(EVAL(3)/scalrd,EVAL(1)/scalrd,EVAL(2)/scalrd) 
EL3=rscalrd*RD(EVAL(1)/scalrd,EVAL(2)/scalrd,EVAL(3)/scalrd) 

DO 70 1=1,6 
DO 70 J=l,6 

70 RF(I,J)=0. 

RF (2,1)=KSC*EL1 
RF (4,3)=KSC*EL2 
RF (6,5)=KSC*EL3*gsq 

C *********************************************************************************** 

c no rotation required 

60 IF (STRATE)THEN 

DO 40 1=1,6 
DO 40 J=l,6 

40 FT (I, J) =F (I, J) +RF (I, J) 
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ELSE 


■ k k k k k k k 


0 

c rotate to lab frame 

c kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk'. 
DO 15 1=1,6 
DO 15 J=l,6 

TEMP(I,J)=0. 

DO 15 M=l,6 

TEMP(I,J)=TEMP(I,J)+RF(I,M)*ROT6(J,M) 
DO 20 1=1,6 
DO 20 J=l,6 

FT (I, J) =F (I, J) 

DO 20 M=1, 6 

FT(I,J)=FT(I,J)+ROT6(I,M)*TEMP(M,J) 
END IF 
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20 


■ k k k k k k 


'■J l C'> l C'> l C'J l C'j'r'J l C'J l C'J l C'J l C'J l C'jk''J l C'> l C'> l C'5 l C'> l C'5 l C'> l C'> l C 


■■j | c'j | c'j | c'j'r'j | c'j'r'j | c'j | c'> | c'> | c'> | c'5 | c'? | c'> | c'> | c'5 | c'> | c'j | c'j'r'j | c'j | c'j | c'j | c'j | c'jk''j | c'? | c'> | c'> | c'? | c'? | c'> | c'> | c 


C icic'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k: 

c END SPACE CHARGE SECTION 
c **************************: 
c START NO SPACE CHARGE SECTION 

ELSE lie. KSC=0. 

DO 23 1=1,6 
DO 23 J=l,6 

23 FT (I, J) =F (I, J) 

END IF 

c END NO SPACE CHARGE SECTION 


c * * * * - 


it'k'k'kir'k'k'k'k'k' 


' -k -k -k -k -k -k ' 


k'kk'kk'kk'. 


C MATRIX OPERATIONS TO OBTAIN DIFF. EQUATIONS 


c kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
DO 30 1=1,6 
DO 30 J=l,I 
DSX(I,J)=0. 

DO 25 M=1,6 

2 5 DSX (I, J) =DSX (I, J) +FT (I, M) *SX (M, J) +SX (I, M) *FT (J,M) 

30 DSX(J,I)=DSX(I,J) 

DO 26 1=7,12 
DO 26 J=l,6 
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DSX(I,J)=0. 

DO 2 6 M=l,6 

26 DSX(I,J)=DSX(I,J)+FT(1-6,M)*SX(M+6,J) 

RETURN 
END 


4 Example 


The TRIUMF injector electron linac, EINJ, takes bunches from energy of 
300 keV to ~ 10 MeV if properly phased and the peak gradient is 20 MV/m. 
Below is example for phase 6 — 0 at the start of the calculation. 



(c) R. Baartman 05/26/15 distance (cm) 

Red is the 2rms transverse size, and green is the 2rms longitudinal (bunch 
length). The input bunch parameters are somewhat arbitrary, roughly the 
condition for a minimum beam size at exit. This particular case has zero bunch 
charge. 


15 







































In this second example, TRANSOPTR is instructed to fit the 65 matrix element to 
zero. This makes energy insensitive to input phase, thus finding the peak energy 
gain phase. This phase turns out to be 6 = —15.46°. 



In the third example, bunch charge has been raised to 30 pC. 
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4.1 Timing 

Each calculation above takes roughly 400 Runge-Kutta steps for 2400 calls to the 
SCLINAC routine. This gives 5-figure accuracy to the transfer matrix and the 
a-matrix, and is easily enough for describing reality considering that the on-axis 
field is only known to 2 or 3 significant figures. The extra accuracy is useful, 
however for fitting matrix or beam matching, which is done with a downhill 
simplex method, or simulated annealing for cases of more than 3 fitting 
parameters. 

On my unremarkable, circa 2006 Intel desktop, each run through the linac takes 
about 17 milliseconds with zero bunch charge and 25 milliseconds with space 
charge. The difference is due to the Carlson elliptic integrals needed for the 
space charge case. 

On a typical optics matching case, one varies 2 solenoids, the buncher amplitude, 
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and the linac phase, to minimize the bunch size and energy spread at the linac 
output. A calculation with such a fit requires typically a half million total calls to 
SC (the space charge routine for no-linac case) and SCLINAC, and so takes about 
5 seconds CPU time. The result is shown below. The bunch charge is 15 pC. 

Rather extravagantly, each calculation starts from the cathode whereas it would 
have been more efficient to store the beam parameter set at the buncher entrance 
and start it from there. The DC acceleration to 300 keV from the cathode is 
described in reference 0. 

The Buncher itself, located at s = 85 cm, is calculated as just another linac, 
phased to give no energy gain. 
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The code describing the whole transport from cathode to linac exit, is shown 
below. 


SUBROUTINE TSYSTEM 

COMMON /BLOC1/cO,Ebun,phbun,cl,c2,phase,grad,birk 
COMMON/MOM/P,BRHO,pMASS,ENERGK,GSQ,ENERGKi,charge,current 
COMMON/SS/SX(13,6) 
nss=2 
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c sollen=6.32 this was with old clamps 

sollen=5.64 
sollenrot = 7.2 9 
solt=c0*140.0/10000 . 
soll=cl*116.9/10000 . 
sol2=c2*116.9/10000 . 
c following are from autocad dwg. 
atgunsol=42.193 
atdb=5 6.9 
atrfgap=85.642 
atdb2=105.7 
atsoll = 122.777 
atsol2 = 17 6.577 
atmatch=237. 
sap=2.5/2.*2.54 
sss=sollen/nss 

c since actually a solenoid's aberration is unrelated to its length, 
wt=ww/sqrt(float(nss) ) 
acclen=12.0 

call axez (59, 121,-100.,acclen, 10, 1,2) 
call drs(atgunsol-acclen-sollen/2.) 
rho=100.*brho/solt 
th= sss/rho/2.*57.29578 
do i=l,nss 

if(i.eq.(nss+l)/2)then 
call solenoidn(solt,sss,sap,wt,'SOLg') 
else 

call solenoidn(solt,sss,sap,wt,' .') 

endif 

call rotate(-th,' .') 

enddo 

th= sollenrot/rho/2.*57.29578 
call rotate(th,' .') 

if(iprint.ne.0)write(6,*)'Sol.B/1gauss,rot.= ',solt*10000.,th*ns 
call drs(atdb-atgunsol-sollen/2.) 
call drs(atrfgap-atdb-10.) !-ccy 

call linacn(62,21,Ebun,20.,1.3e9,phbun,20,'BNCH') iBuncherEH 
sollen=5.85 
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sollenrot=7.37 

sss=sollen/nss 

call drs(atdb2-atrfgap-10.) 

call drs(atsoll-atdb2-sollen/2.) 

rho=100.*brho/soli 

th= sss/rho/2.*57.29578 

do i=l,nss 

if(i.eq.(nss+1)/2)then 
call solenoidn (soil,sss,sap, wt, 'SOLI'' ) 
else 

call solenoidn(sol1,sss,sap,wt, ' .') 

endif 

call rotate(-th ,' .') 

enddo 

th= sollenrot/rho/2.*57.29578 
call rotate(th ,' .') 

if(iprint.ne.O)write(6 ,*)' Sol.B/1gauss,rot.= ',3011*10000.,th*ns 

call drs(atsol2-atsoll-sollen) 

rho=100.*brho/sol2 

th= sss/rho/2.*57.2957 8 

do i=l,nss 

if(i.eq.(nss+1)/2)then 
call solenoidn(sol2,sss,sap,wt,'SOL2') 
else 

call solenoidn(sol2,sss,sap,wt, ' .') 

endif 

call rotate(-th,' .') 

enddo 

th= sollenrot/rho/2.*57.29578 
call rotate(th,' .') 

if(iprint.ne.O)write(6,*)'Sol.B/1gauss,rot.= ',sol2*10000.,th*ns 
d5=74.96-sollen/2. 

call drs(atmatch-atsol2-sollen/2.) 

call drift (0.,'Mate') 

call drs (15.-1.023) 

call fit (1,1, 1,0 ., 1., 1) 

call fit(1,5,5,0.4,8.,1) 

call drift (1.023,' . ' ) 
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sx(13,5)=0. 

call linacn(71,447,grad,127.953,1.3e9,phase,100,'EINJ') 

call drift(1.023, ' .' ) 

call fit (1,1,1,0 ., 4 ., 1) 

call fit(1,2,2,0., 1000 . , 1) 

call fit(1,5,5,0.,70.,1) 

call fit (1, 6, 6, 0., 1000 . , 1) 

enrg=sx(13,6)*0.511 

c we don't want output energy to be too far off. 
call fitarb(10.02,enrg,100.,1) 
call drs(61.37) 
call drift(0.,'MdbO') 
call vective (1) 
return 
end 
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